(CC-BY-NC-SA)
This document was produced on Mo Sep 28 2015 using mapview version 0.7.3
The easiest way to install mapview is to use library("devtools"). With devtools installed and loaded type the following:
For the stable release of mapview:
install_github("environmentalinformatics-marburg/mapview")
For the devolpment version:
install_github("environmentalinformatics-marburg/mapview", ref = "develop")
Working with spatial data in R I find myself quite often in the need to quickly visually check whether a certain analysis has produced reasonable results. There are two ways I usually do this. Either I:
Both these approaches are semi-optimal. Where option 1. is fine for a quick glance at a coarse patterns, it lacks the possibility to have a closer look into the results via zooming and panning. While option 2. provides the interactivity, the detour via the hard disk is annoying (at best), especially when fine-tuning and checking regularly.
I attended this years useR2015! conference in Aalborg (which was marvelous!) and attended the session on interactive graphics in R where Joe Cheng from RStudio presented the leaflet package. leaflet is great as it gives you a great deal of control over the individual map components. This, however, also means that in order to get some useful visualization of spatial objects, we need to do a fair bit of coding for the map components to show what we want (e.g. we would need to provide popup-text manually for each object that we want to map). What a GIS-like functionality would need is some default behaviour for different objects from the spatial universe.
This got me thinking and sparked my enthusiasm to write some wrapper functions for leaflet to provide at least very basic GIS-like interactive graphing capabilities that are directly accessible within RStudio (or the web browser, if you’re not using RStudio).
The result is a package called mapview. The main workhorse function is mapView() (mind the capital ‘V’) and is currently defined for:
A call to mapView() will return an object of class mapview. This class has 2 slots
leaflet map. This is an S3 class object (see leaflet documentation for details on the specifics).By default mapView() provides two base layers between which one can toggle (a standard OpenStreeMap layer and ESRI’s WorldImagery layer). Depending on the object and/or argument settings one or several layers are created (each with or without it’s own legend) that can also be toggled. Note that in order to render properly, all layers need to be re-projected to leaflet’s underlying web mercator projection
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
which in the case of large Raster* objects or Spatial objects with lots of features can be time consuming.
In the following I would like to present a few use case scenarios that highlight the current capabilities of mapview (largely taken from the help files):
NOTE: similar to spplot() for Raster* objects mapView() has a maxpixels = argument to avoid long rendering times for large rasters. This is currently set to a default value of 500000 (which produces acceptable rendering times on my machine) and a suitable warning is printed for larger Raster objects.
library(mapview)
library(raster)
data(meuse.grid)
coordinates(meuse.grid) = ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) = TRUE
## convert to RasterStack
meuse_rst <- stack(meuse.grid)
mapView(meuse_rst[[4]])
Actually, in the original meuse.grid data frame the 4th column is a factor with 3 levels. So when we convert the RasterLayer to a factor, we will see that we get a proper factorial layer with a suitable legend.
meuse_rst[[4]] <- as.factor(meuse_rst[[4]])
mapView(meuse_rst[[4]])
If you pass a RasterStack/Brick to mapView() this will create one map layer + legend for each RasterLayer. This is likely to cause problems for larger Stacks as there is currently no way to tell mapView() to only show the legend of the highlighted layer. This means that some of the legends won’t show in the viewer.
mapView(meuse_rst)
There is (to my knowledge) currently no way to provide this functionality in the leaflet package, though I have seen solutions for this in JavaScript which means that this will likely be available in the future.
Given that a SpatialPixelsDataFrame has an attribute table we can use argument zcol = "column_name" to plot specific columns of the table. The zcol argument can be used with the …DataFrame versions of all Spatial* objects that are currently supported (see below).
mapView(meuse.grid, zcol = "dist")
For spatial vector objects from the sp-package I would like to highlight the following arguments:
zcol similar to SpatialPixelsDataFrame you can specify which column of the attribute table to view. If this is set to NULL (the default) all columns will be used.burst if set to FALSE (the default) one layer is rendered with information from all columns in the attribute table provided in pop-ups (triggered by clicking on the features). If set to TRUE one layer + legend for each column in the attribute table will be drawn. Again, depending on the number of columns, this may mean that not all legend will be visible.data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
# only one layer, all info in popups
mapView(meuse)
Another argument you can use with SpatialPointsDataFrame is radius which takes either an integer value (default is 10) or the character name of one of the columns in the attribute table. If the latter is provided, circles are scaled relative to the respective values in the column provided.
mapView(meuse, radius = "lead")
Apart from the radius argument, the same functionality is implemented for
data("gadmCHE")
mapView(gadmCHE)
and
data("atlStorms2005")
mapView(atlStorms2005, burst = TRUE)
For the remote sensing inclined, mapview-package offers some additional functionality:
viewRGB() to view true-/false color images of satellite imageryviewExtent() to only view the extent of a Raster* object. This is useful to quickly check the area covered by an image. This is really fast even for large rasters, as only the extent (4 coordinates) need to be re-projected. Note, viewExtent() also works for all supported Spatial objects.viewRGB() provides the possibility to view Red-Green-Blue true-/false-color images of RasterStacks/Bricks. Any band combination is possible.
data("poppendorf")
viewRGB(poppendorf, 4, 3, 2) #Bands 4,3,2 correspond to Red Green Blue
Note, by default the color to render NA values is set to "transparent", so let’s change that to "black"
viewRGB(poppendorf, 5, 4, 3, na.color = "black") #Bands 5,4,3 correspond to NIR Red Green
Fun fact: within this tiny area one can find at least five breweries (see below)… Welcome to the land of Beer!
As mentioned above, it is possible to view only the extent of Raster and bbox of Spatial objects. Clicking on the rectangle triggers a pop-up window showing the coordinates of the corners of the extent in epsg:4326.
viewExtent(poppendorf)
viewExtent(meuse, color = "black")
To make mapview even more user-friendly, I have defined a +-method for combinations of
mapview + mapviewmapview + ANY (ANY here refers to all classes supported by mapView())leaflet + ANY (so that one can easily pass spatial objects to existing leaflet maps)This means you can easily combine objects by enabling things like
mapView(meuse.grid, zcol = "soil") + meuse
Basically, any combination is possible.
The zoom will be set to the extent of the objected that was added last. Also, layers are overlaid according to the order of calls. Note, however, that leaflet will always render vector data on top of raster data,
data("breweries91")
mapView(breweries91, color = "red") + viewRGB(poppendorf, 4, 3, 2) + viewExtent(poppendorf)
So here you have it, 5 breweries within the extent of the RGB image
It is possible to view objects that are not projected (i.e. have their coordinate reference system - CRS - set to NA). If this is the case, coordinates are rescaled to lie between 0 and 1 and then the y-coordinates are multiplied by the ratio of y-extent / x-extent. This means that x-coordinates will always have min == 0 and max == 1 while y-coordinates will have min == 0 but their max will depend on the ratio, so that this can be slightly higher or lower than 1.
When we supply non-projected data to mapView() a suitable warning will be printed to the console (see below).
Here’s an example for point data. For this we load the meuse data set once again, but this time we won’t set the CRS.
data("meuse")
coordinates(meuse) <- ~x+y
## check coordinate reference system
proj4string(meuse)
## [1] NA
Yet, even though it has no CRS, we can still visulize it. Though we obviously cannot provide any background maps for the visualization.
mapView(meuse)
## Warning in spCheckAdjustProjection(x): supplied SpatialPointsDataFrame has no projection information!
## scaling coordinates and showing layer without background map
Note, how the coordinates have been rescaled.
This is also possible for Raster* objects, which allows us to view pretty much any type of raster data, including those that cannot be projected, such as photos. This is quite handy as it enables closer investigation of non-standard spatial data, for example gathered through imaging ground-based remote sensing techniques (e.g. thermal imaging).
Here’s a toy example using one of the header images from the Environmental Informtics Marburg web page that shows the author of this package during an R teaching session at our research station at Mt. Kilimanjaro in Tanzania.
## get image from the web page
library(jpeg)
web_image <- "http://umweltinformatik-marburg.de/uploads/tx_rzslider/teaching_header_kili_resized.jpg"
jpg <- readJPEG(readBin(web_image, "raw", 1e6))
## Convert image data to raster
rst_blue <- raster(jpg[, , 1])
rst_green <- raster(jpg[, , 2])
rst_red <- raster(jpg[, , 3])
img <- brick(rst_red, rst_green, rst_blue)
## visualize image in viewer
viewRGB(img)
## Warning in rasterCheckAdjustProjection(x, maxpixels): supplied RasterBrick has no projection information!
## scaling coordinates and showing layer without background map
In case you are interested, under the hood, the rendering of unprojected Raster* data works similar to that of Spatial* objects. Corner coordinates of the extent are rescaled (y-coordinates are multiplied with the ratio of nrow / ncol to preserve dimensions). In this case, however, we need to provide a CRS as leaflet expects one. Therefore, we set this to standard latitude/longitude EPSG:4236. This means, that the rendered image actually does have a spatial reference but mapview simply does not provide basemaps. If we were to provide one, we would see that our data lies somewhere in the Gulf of Guinea.
Applications for non-projected data are obviously limited, yet, especially the possibility to view standard images interactively provides a nice tool for checking parts of an image(stack) in more detail.
Keep in mind that speed and memory limitations for Raster* data still applies, though speed will be much less of an issue, as we do not need to reporject the data.
Integrating mapview and leaflet works both ways…
You can first create a map using leaflet and then add spatial objects as you like with mapview allowing you to flexibly create maps in whatever way you like.
## first set the CRS again
proj4string(meuse) <- CRS("+init=epsg:28992")
m <- leaflet() %>%
addProviderTiles("Stamen.Toner")
m + meuse
Given that objects created with mapView() contain a slot called map, you can first create a map using mapView() and then add components provided by leaflet. To reproduce the following example you will (at the time of this writing) need to install leaflet from the GitHub leaflet repository which provides new functionality to measure distances and areas. The following code will then add a control element in the top right corner of the map which lets you interactively create a polygon and will show its perimeter and area. Double-click to finish a measurement.
m <- mapView(breweries91) + viewExtent(poppendorf)
m@map %>% addMeasure(primaryLengthUnit = "meters", primaryAreaUnit = "sqmeters",
activeColor = "#3D535D", completedColor = "#7D4479")
So now you can measure the distances between the breweries or figure out how big the area is that is covered by the extent…
As a final, more advanced, example of what can be done with mapview we will mimic a common workflow when working with spatial data. We will calculate the mean January temperature for each district of Switzerland using raster::extract. Then we will render this on top of the mean January temperatures and also add a hill-shade raster of Switzerland. All data sets are included in mapview. They were created following the examples shown in the help pages of raster::getData and raster::hillShade. This example also highlights a few arguments that can be adjusted which were not mentioned so far.
data("tmeanCHE")
data("hillshadeCHE")
data("gadmCHE")
t <- extract(tmeanCHE, gadmCHE, fun = mean, na.rm = TRUE, sp = TRUE)
t$Jan.mean.T <- round(t$Jan.mean.T, 2)
m1 <- mapView(hillshadeCHE, color = grey.colors(3),
layer.opacity = 1, layer.name = "Hillshade",
map.types = "Acetate.hillshading", legend = FALSE)
m2 <- mapView(tmeanCHE, layer.opacity = 0.6,
map.types = c("OpenStreetMap.BlackAndWhite",
"OpenStreetMap.HOT"),
layer.name = "JanTemp")
m3 <- mapView(t, zcol = "Jan.mean.T", fillOpacity = 0.9)
m1 + m2 + m3
Note how we
fillOpacity to the underlying leaflet call for the polygons layerI hope this highlights that despite providing sensible default behaviour for viewing spatial data mapview also provides a fair deal of flexibility.
… still wonder where exactly the difference between mapview and leaflet lies, consider the following:
Let’s reproduce a simple mapView() call using leaflet()
Here’s the mapview version:
mapView(meuse)
And now let’s get the identical result using leaflet:
## first we need to reproject meuse to geographic coordinates
meuse <- spTransform(meuse,
CRSobj = CRS("+proj=longlat +datum=WGS84 +no_defs"))
## then we need to generate the text for the pop-ups
df <- as.data.frame(sapply(meuse@data, as.character),
stringsAsFactors = FALSE)
nms <- names(df)
txt_x <- paste0("x: ", round(coordinates(meuse)[, 1], 2))
txt_y <- paste0("y: ", round(coordinates(meuse)[, 2], 2))
txt <- rbind(sapply(seq(nrow(meuse@data)), function(i) {
paste(nms, df[i, ], sep = ": ")
}), txt_x, txt_y)
txt <- sapply(seq(ncol(txt)), function(j) {
paste(txt[, j], collapse = " <br/> ")
})
## finally we can create our map
leaflet(meuse) %>%
addProviderTiles("OpenStreetMap", group = "OpenStreetMap") %>%
addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery") %>%
addCircleMarkers(lng = coordinates(meuse)[, 1],
lat = coordinates(meuse)[, 2],
group = "meuse",
color = mapViewPalette(3)[length(mapViewPalette(3))],
popup = txt) %>%
addLayersControl(position = "bottomleft",
baseGroups = c("OpenStreetMap",
"Esri.WorldImagery"),
overlayGroups = "meuse")
I guess this highlights quite well what mapview is intended for. It serves to save you some effort to interactively examine your spatial data by giving you the possibility to render it using a one-liner of code.
Where to go from here?
Well, there are still a lot of limitations that need addressing:
In future releases I would like to
I hope that mapview will prove useful for some people out there.
If you have any feedback, please don’t hesitate to contact me.
Bug reports should be filed at https://github.com/environmentalinformatics-marburg/mapview/issues
Best,
Tim